function [L,U] = gauss(A)

% ELIMINAZIONE DI GAUSS
% Dati di input:
%   A = matrice quadrata
% Dati di output:
%   L = matrice triangolare inferiore
%   U = matrice triangolare superiore
% tali che A = LU

[n,m] = size(A);
if n ~= m
   disp('ERRORE: la matrice di input non e'' quadrata')
   break
end

for k=1:n-1

% controllo sulla grandezza dell'elemento diagonale:
   if abs(A(k,k)) <  realmin
      disp('Un elemento diagonale e'' piccolo')
      warning('La fattorizzazione potrebbe essere inaccurata')
   end

   A(k+1:n,k) = A(k+1:n,k)/A(k,k);
   A(k+1:n,k+1:n) = A(k+1:n,k+1:n) - ...
   A(k+1:n,k)*A(k,k+1:n);

end

L = tril(A,-1) + eye(n);
U = triu(A);
